Setup

Clear previous, load packages, manage error

load and prep variables

reference for clusters and summary of mixed effects models for each cluster

-Considered running ROI as a factor, but then need a reference cluster.
-It also seems like a reasonable assumption that the clusters would show different functional forms over time. Task activation clusters.

course description with box plots

-no extreme outliers when averaging across sessions
-frequent and infrequent are highly similar

ggplot(clpsc, aes(x=roi,y=psc,fill=cond)) + geom_violin(alpha=0.7,draw_quantiles=c(0.5)) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + facet_wrap(~group)

ggplot(clpsc, aes(x=roi,y=psc,fill=cond)) + geom_boxplot(outlier.colour = "#1F3552", outlier.shape = 20) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + facet_wrap(~group)

ggplot(clpsc, aes(x=roi,y=psc,fill=group)) + geom_boxplot(outlier.colour = "#1F3552", outlier.shape = 20) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + facet_wrap(~cond)

plot individuals over session for each roi

-frequent and infrequent highly similar
-linear model with quadratic term would work all rois

#zfclus1_lifg
g1 <- ggplot(data = subset(clpsc,roi=="zfclus1_lifg"), aes(x = run, y = psc, colour=cond))
g2 <- g1 + geom_line() + geom_point() + facet_wrap(~sub)
g3 <- g2 + theme_bw() + scale_x_continuous(name = "run")
g4 <- g3 + scale_y_continuous(name = "psc") + ggtitle("zfclus1_lifg")
print(g4)

#zfclus1_lins
g1 <- ggplot(data = subset(clpsc,roi=="zfclus1_lins"), aes(x = run, y = psc, colour=cond))
g2 <- g1 + geom_line() + geom_point() + facet_wrap(~sub)
g3 <- g2 + theme_bw() + scale_x_continuous(name = "run")
g4 <- g3 + scale_y_continuous(name = "psc") + ggtitle("zfclus1_lins")
print(g4)

#zfclus2
g1 <- ggplot(data = subset(clpsc,roi=="zfclus2"), aes(x = run, y = psc, colour=cond))
g2 <- g1 + geom_line() + geom_point() + facet_wrap(~sub)
g3 <- g2 + theme_bw() + scale_x_continuous(name = "run")
g4 <- g3 + scale_y_continuous(name = "psc") + ggtitle("zfclus2")
print(g4)

#zfclus3
g1 <- ggplot(data = subset(clpsc,roi=="zfclus3"), aes(x = run, y = psc, colour=cond))
g2 <- g1 + geom_line() + geom_point() + facet_wrap(~sub)
g3 <- g2 + theme_bw() + scale_x_continuous(name = "run")
g4 <- g3 + scale_y_continuous(name = "psc") + ggtitle("zfclus3")
print(g4)

#zfclus4_brainstem
g1 <- ggplot(data = subset(clpsc,roi=="zfclus4_brainstem"), aes(x = run, y = psc, colour=cond))
g2 <- g1 + geom_line() + geom_point() + facet_wrap(~sub)
g3 <- g2 + theme_bw() + scale_x_continuous(name = "run")
g4 <- g3 + scale_y_continuous(name = "psc") + ggtitle("zfclus4_brainstem")
print(g4)

#zfclus4_brainstem_aan_dr
g1 <- ggplot(data = subset(clpsc,roi=="zfclus4_brainstem_aan_dr"), aes(x = run, y = psc, colour=cond))
g2 <- g1 + geom_line() + geom_point() + facet_wrap(~sub)
g3 <- g2 + theme_bw() + scale_x_continuous(name = "run")
g4 <- g3 + scale_y_continuous(name = "psc") + ggtitle("zfclus4_brainstem_aan_dr")
print(g4)

#zfclus4_brainstem_aan_ppn
g1 <- ggplot(data = subset(clpsc,roi=="zfclus4_brainstem_aan_ppn"), aes(x = run, y = psc, colour=cond))
g2 <- g1 + geom_line() + geom_point() + facet_wrap(~sub)
g3 <- g2 + theme_bw() + scale_x_continuous(name = "run")
g4 <- g3 + scale_y_continuous(name = "psc") + ggtitle("zfclus4_brainstem_aan_ppn")
print(g4)

#zfclus4_parahipp
g1 <- ggplot(data = subset(clpsc,roi=="zfclus4_parahipp"), aes(x = run, y = psc, colour=cond))
g2 <- g1 + geom_line() + geom_point() + facet_wrap(~sub)
g3 <- g2 + theme_bw() + scale_x_continuous(name = "run")
g4 <- g3 + scale_y_continuous(name = "psc") + ggtitle("zfclus4_parahipp")
print(g4)

#zfclus5
g1 <- ggplot(data = subset(clpsc,roi=="zfclus5"), aes(x = run, y = psc, colour=cond))
g2 <- g1 + geom_line() + geom_point() + facet_wrap(~sub)
g3 <- g2 + theme_bw() + scale_x_continuous(name = "run")
g4 <- g3 + scale_y_continuous(name = "psc") + ggtitle("zfclus5")
print(g4)

#zfclus6
g1 <- ggplot(data = subset(clpsc,roi=="zfclus6"), aes(x = run, y = psc, colour=cond))
g2 <- g1 + geom_line() + geom_point() + facet_wrap(~sub)
g3 <- g2 + theme_bw() + scale_x_continuous(name = "run")
g4 <- g3 + scale_y_continuous(name = "psc") + ggtitle("zfclus6")
print(g4)

prep variables for modeling

-center run like in behavioral analyses
-when centering time on run, linear effects of run reflect slope at intercept which is at the middle of the day

clpsc$run0 <- clpsc$run-1
clpsc$runc3 <- clpsc$run-3
clpsc$sub <- as.factor(clpsc$sub)
clpsc$group <- as.factor(clpsc$group)
clpsc$roi <- as.factor(clpsc$roi)
clpsc$cond <- as.factor(clpsc$cond)

set up contrasts for group and cond

levels(clpsc$group)
## [1] "older"   "younger"
options(contrasts=c("contr.sum", "contr.poly"))
contrasts(clpsc$group)  #older = +1; younger = -1
##         [,1]
## older      1
## younger   -1
contrasts(clpsc$cond) #freq = +1; infreq = -1
##        [,1]
## freq      1
## infreq   -1
#nothing in the data supports that Infrequent is any different so exclude from model for simplicity to match with activation analyses

for each model, start with most complex and trim back

-more complex models are commented out for simplicity
-To comment and uncomment blocks of code, use ctrl+shift+c

zfclus1_lifg (LIFG)

# #full
# model1 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus1_lifg"))
# summary(model1)
# 
# model1b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus1_lifg"))
# summary(model1b)


# #drop random quadratic
# model2 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus1_lifg"))
# summary(model2)
# 
# model2b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus1_lifg"))
# summary(model2b)
# 
# anova(model1,model2)
# #non-sig, drop rand quadratic
# #
# 
# anova (model1b,model2b)
# #non-sig, drop rand quadratic

# #drop random slopes
# model3 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 | sub), data = subset(clpsc,roi=="zfclus1_lifg"))
# summary(model3)
# anova(model2,model3)
# #non-sig, drop random slopes
# 
# model3b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 | sub), data = subset(clpsc,roi=="zfclus1_lifg"))
# summary(model3b)
# anova(model2b,model3b)
# #non-sig, drop random slopes


# #remove fixed quadratic
# model4 <- lmer(psc ~ runc3*group*cond + (1 | sub), data = subset(clpsc,roi=="zfclus1_lifg"))
# summary(model4)
# anova(model3,model4)
# #non-sig, keep only random intercepts

model4b <- lmer(psc ~ runc3*group + (1 | sub), data = subset(clpsc,roi=="zfclus1_lifg"))
summary(model4b)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: psc ~ runc3 * group + (1 | sub)
##    Data: subset(clpsc, roi == "zfclus1_lifg")
## 
## REML criterion at convergence: -94.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5315 -0.5811  0.0044  0.6018  3.3766 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  sub      (Intercept) 0.02802  0.1674  
##  Residual             0.02849  0.1688  
## Number of obs: 240, groups:  sub, 24
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)    0.002763   0.035863  22.000000   0.077 0.939281    
## runc3         -0.037630   0.007704 214.000000  -4.884 2.03e-06 ***
## group1         0.055416   0.035863  22.000000   1.545 0.136558    
## runc3:group1  -0.030504   0.007704 214.000000  -3.959 0.000102 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) runc3 group1
## runc3       0.000              
## group1      0.000  0.000       
## runc3:grop1 0.000  0.000 0.000
# anova(model3b,model4b)
# #non-sig, keep only random intercepts

# anova(model4a,model4b)

final_model=model4b

#plot by groups and day
g1 <- ggplot(data = subset(clpsc,roi=="zfclus1_lifg"), aes(x = runc3, y = psc, shape = group, color=group))
g2 <- g1 + stat_summary(fun.data=mean_se, geom="pointrange", size=.5) + scale_colour_manual(values=c("#6666FF","#666666"))
g3 <- g2 + stat_summary(aes(y = fitted(final_model), linetype = group,color=group), fun.y = mean, geom="line")
g4 <- g3 + theme_bw(base_size=18) 
g5 <- g4 + scale_y_continuous(name = "Percent Signal Change") + scale_x_continuous(label=c("1","2","3","4","5"), name = "Session")
print(g5)

ggsave(filename="zfclus1_lifg.png")
## Saving 7 x 5 in image
#example to write out raw and fitted (see rachel's latest code)
#writeout=cbind(fitted(final_model),subset(clpsc,roi=="zfclus1_lifg"))
#write.csv(writeout,"final_model_zfclus1_lifg.csv",row.names=FALSE, na="")

can output individual estimates of slope via random effects

-but careful because n=12 per group is likely not enough for reliable individual differences analyses

model4a <- lmer(psc ~ runc3*group + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus1_lifg"))
summary(model4a)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: psc ~ runc3 * group + (1 + runc3 | sub)
##    Data: subset(clpsc, roi == "zfclus1_lifg")
## 
## REML criterion at convergence: -97.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5771 -0.5679 -0.0153  0.5713  3.1233 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  sub      (Intercept) 0.0281467 0.16777       
##           runc3       0.0006236 0.02497  -0.56
##  Residual             0.0272098 0.16495       
## Number of obs: 240, groups:  sub, 24
## 
## Fixed effects:
##               Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   0.002763   0.035863 22.000000   0.077  0.93928    
## runc3        -0.037630   0.009092 22.000000  -4.139  0.00043 ***
## group1        0.055416   0.035863 22.000000   1.545  0.13656    
## runc3:group1 -0.030504   0.009092 22.000000  -3.355  0.00286 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) runc3  group1
## runc3       -0.298              
## group1       0.000  0.000       
## runc3:grop1  0.000  0.000 -0.298
coef(model4a)
## $sub
##      (Intercept)        runc3     group1 runc3:group1
## 100 -0.093332782 -0.029548362 0.05541605  -0.03050375
## 102 -0.150789161 -0.026969171 0.05541605  -0.03050375
## 106  0.202334730 -0.058754884 0.05541605  -0.03050375
## 108  0.086445024 -0.044051740 0.05541605  -0.03050375
## 112 -0.280631843 -0.024436509 0.05541605  -0.03050375
## 113  0.031207244 -0.055589624 0.05541605  -0.03050375
## 115 -0.145313544 -0.023921959 0.05541605  -0.03050375
## 117  0.162352724 -0.060079342 0.05541605  -0.03050375
## 118 -0.238431632 -0.007232585 0.05541605  -0.03050375
## 119  0.049107089 -0.051085053 0.05541605  -0.03050375
## 120  0.091195921 -0.018108192 0.05541605  -0.03050375
## 121  0.319014630 -0.051784903 0.05541605  -0.03050375
## 200  0.143512951 -0.048426689 0.05541605  -0.03050375
## 203 -0.120642515 -0.023439342 0.05541605  -0.03050375
## 204  0.227700372 -0.062436403 0.05541605  -0.03050375
## 205  0.127581181 -0.051903250 0.05541605  -0.03050375
## 207  0.067226579 -0.044869232 0.05541605  -0.03050375
## 208  0.017615916 -0.025633391 0.05541605  -0.03050375
## 209 -0.006369182 -0.032065803 0.05541605  -0.03050375
## 210  0.077653366 -0.056301472 0.05541605  -0.03050375
## 211 -0.219329425 -0.010221398 0.05541605  -0.03050375
## 214 -0.102815828 -0.020174621 0.05541605  -0.03050375
## 216 -0.083929089 -0.036414756 0.05541605  -0.03050375
## 217 -0.095045927 -0.039675967 0.05541605  -0.03050375
## 
## attr(,"class")
## [1] "coef.mer"

zfclus1_lins (LINS)

# #full
# model1 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus1_lins"))
# summary(model1)
# 
# model1b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus1_lins"))
# summary(model1b)

#drop random quadratic
# model2 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus1_lins"))
# summary(model2)
# 
# model2b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus1_lins"))
# summary(model2b)
# 
# anova(model1,model2)
# #non-sig, drop rand quadratic
# 
# anova(model1b,model2b)
# #non-sig, drop rand quadratic

# #drop random slopes
# model3 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 | sub), data = subset(clpsc,roi=="zfclus1_lins"))
# summary(model3)
# anova(model2,model3)
# #non-sig, drop random slopes

model3b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 | sub), data = subset(clpsc,roi=="zfclus1_lins"))
summary(model3b)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: psc ~ runc3 * group + I(runc3^2) * group + (1 | sub)
##    Data: subset(clpsc, roi == "zfclus1_lins")
## 
## REML criterion at convergence: -48.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6694 -0.5471 -0.0805  0.5544  3.6294 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  sub      (Intercept) 0.03587  0.1894  
##  Residual             0.03219  0.1794  
## Number of obs: 240, groups:  sub, 24
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)         0.244809   0.042665  27.440000   5.738 3.99e-06 ***
## runc3              -0.057674   0.008190 212.000000  -7.042 2.59e-11 ***
## group1              0.079045   0.042665  27.440000   1.853  0.07471 .  
## I(runc3^2)          0.019482   0.006921 212.000000   2.815  0.00534 ** 
## runc3:group1       -0.023530   0.008190 212.000000  -2.873  0.00448 ** 
## group1:I(runc3^2)  -0.013178   0.006921 212.000000  -1.904  0.05828 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) runc3  group1 I(3^2) rnc3:1
## runc3        0.000                            
## group1       0.000  0.000                     
## I(runc3^2)  -0.324  0.000  0.000              
## runc3:grop1  0.000  0.000  0.000  0.000       
## grp1:I(3^2)  0.000  0.000 -0.324  0.000  0.000
# anova(model2b,model3b)
# #non-sig, drop random slopes


#remove fixed quadratic
# model4 <- lmer(psc ~ runc3*group*cond + (1 | sub), data = subset(clpsc,roi=="zfclus1_lins"))
# summary(model4)
# 
# model4b <- lmer(psc ~ runc3*group + (1 | sub), data = subset(clpsc,roi=="zfclus1_lins"))
# summary(model4b)
# 
# anova(model3,model4)
# #sig, keep fixed quadratic
# 
# anova(model3b,model4b)
# #sig, keep fixed quadratic


final_model=model3b

#plot by groups and day
g1 <- ggplot(data = subset(clpsc,roi=="zfclus1_lins"), aes(x = runc3, y = psc, shape = group, color=group))
g2 <- g1 + stat_summary(fun.data=mean_se, geom="pointrange", size=.5) + scale_colour_manual(values=c("#6666FF","#666666"))
g3 <- g2 + stat_summary(aes(y = fitted(final_model), linetype = group,color=group), fun.y = mean, geom="line")
g4 <- g3 + theme_bw(base_size=18) 
g5 <- g4 + scale_y_continuous(name = "Percent Signal Change") + scale_x_continuous(label=c("1","2","3","4","5"), name = "Session")
print(g5)

ggsave(filename="zfclus1_lins.png")
## Saving 7 x 5 in image

zfclus2 (LSLOC)

# #full
# model1 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus2"))
# summary(model1)
# 
# model1b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus2"))
# summary(model1b)

# #drop random quadratic
# model2 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus2"))
# summary(model2)
# 
# model2b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus2"))
# summary(model2b)
# 
# anova(model1,model2)
# #non-sig, drop rand quadratic
# 
# anova(model1b,model2b)
# #non-sig, drop rand quadratic


# #drop random slopes
# model3 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 | sub), data = subset(clpsc,roi=="zfclus2"))
# summary(model3)

# anova(model2,model3)
# #non-sig, drop random slopes

model3b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 | sub), data = subset(clpsc,roi=="zfclus2"))
summary(model3b)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: psc ~ runc3 * group + I(runc3^2) * group + (1 | sub)
##    Data: subset(clpsc, roi == "zfclus2")
## 
## REML criterion at convergence: -98.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2544 -0.5107 -0.0855  0.4737  3.2716 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  sub      (Intercept) 0.04144  0.2036  
##  Residual             0.02518  0.1587  
## Number of obs: 240, groups:  sub, 24
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)         0.323699   0.044514  25.730000   7.272 1.07e-07 ***
## runc3              -0.042415   0.007243 212.000000  -5.856 1.79e-08 ***
## group1              0.030248   0.044514  25.730000   0.680   0.5029    
## I(runc3^2)          0.014043   0.006121 212.000000   2.294   0.0228 *  
## runc3:group1       -0.010562   0.007243 212.000000  -1.458   0.1462    
## group1:I(runc3^2)   0.008912   0.006121 212.000000   1.456   0.1469    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) runc3  group1 I(3^2) rnc3:1
## runc3        0.000                            
## group1       0.000  0.000                     
## I(runc3^2)  -0.275  0.000  0.000              
## runc3:grop1  0.000  0.000  0.000  0.000       
## grp1:I(3^2)  0.000  0.000 -0.275  0.000  0.000
# anova(model2b,model3b)
# # #non-sig, drop random slopes



# #remove fixed quadratic
# model4 <- lmer(psc ~ runc3*group*cond + (1 | sub), data = subset(clpsc,roi=="zfclus2"))
# summary(model4)
# 
# anova(model3,model4)
# #sig, keep fixed quadratic
# 
# 
# model4b <- lmer(psc ~ runc3*group + (1 | sub), data = subset(clpsc,roi=="zfclus2"))
# summary(model4b)
# 
# anova(model3b,model4b)
# #sig, keep fixed quadratic

final_model=model3b

g1 <- ggplot(data = subset(clpsc,roi=="zfclus2"), aes(x = runc3, y = psc, shape = group, color=group))
g2 <- g1 + stat_summary(fun.data=mean_se, geom="pointrange", size=.5) + scale_colour_manual(values=c("#6666FF","#666666"))
g3 <- g2 + stat_summary(aes(y = fitted(final_model), linetype = group,color=group), fun.y = mean, geom="line")
g4 <- g3 + theme_bw(base_size=18) 
g5 <- g4 + scale_y_continuous(name = "Percent Signal Change") + scale_x_continuous(label=c("1","2","3","4","5"), name = "Session")
print(g5)

ggsave(filename="zfclus2.png")
## Saving 7 x 5 in image

zfclus3 (RSLOC)

# #full
# model1 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus3"))
# summary(model1)
# 
# model1b <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus3"))
# summary(model1b)


# #drop random quadratic
# model2 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus3"))
# summary(model2)
# 
# anova(model1,model2)
# #non-sig (marg), drop rand quadratic
# 
# model2b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus3"))
# summary(model2b)
# 
# anova(model1b,model2b)
# #non-sig (marg), drop rand quadratic

# #drop random slopes
# model3 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 | sub), data = subset(clpsc,roi=="zfclus3"))
# summary(model3)
# anova(model2,model3)
# #non-sig, drop random slopes
# 
# model3b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 | sub), data = subset(clpsc,roi=="zfclus3"))
# summary(model3b)
# 
# anova(model2b,model3b)
# #non-sig, drop random slopes
# 
# #remove fixed quadratic
# model4 <- lmer(psc ~ runc3*group*cond + (1 | sub), data = subset(clpsc,roi=="zfclus3"))
# summary(model4)
# anova(model3,model4)
# #non-sig, drop fixed quadratic


model4b <- lmer(psc ~ runc3*group + (1 | sub), data = subset(clpsc,roi=="zfclus3"))
summary(model4b)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: psc ~ runc3 * group + (1 | sub)
##    Data: subset(clpsc, roi == "zfclus3")
## 
## REML criterion at convergence: -68
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.6545 -0.5361 -0.0253  0.5260  2.6920 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  sub      (Intercept) 0.02771  0.1665  
##  Residual             0.03221  0.1795  
## Number of obs: 240, groups:  sub, 24
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)    0.184255   0.035898  22.000000   5.133 3.82e-05 ***
## runc3         -0.057503   0.008192 214.000000  -7.020 2.90e-11 ***
## group1         0.069861   0.035898  22.000000   1.946   0.0645 .  
## runc3:group1  -0.015941   0.008192 214.000000  -1.946   0.0530 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) runc3 group1
## runc3       0.000              
## group1      0.000  0.000       
## runc3:grop1 0.000  0.000 0.000
# anova(model3b,model4b)
# #non-sig, drop fixed quadratic


final_model=model4b

g1 <- ggplot(data = subset(clpsc,roi=="zfclus3"), aes(x = runc3, y = psc, shape = group, color=group))
g2 <- g1 + stat_summary(fun.data=mean_se, geom="pointrange", size=.5) + scale_colour_manual(values=c("#6666FF","#666666"))
g3 <- g2 + stat_summary(aes(y = fitted(final_model), linetype = group,color=group), fun.y = mean, geom="line")
g4 <- g3 + theme_bw(base_size=18) 
g5 <- g4 + scale_y_continuous(name = "Percent Signal Change") + scale_x_continuous(label=c("1","2","3","4","5"), name = "Session")
print(g5)

ggsave(filename="zfclus3.png")
## Saving 7 x 5 in image

zfclus4_brainstem (BSTEM)

# #full
# model1 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus4_brainstem"))
# summary(model1)
# 
# model1b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus4_brainstem"))
# summary(model1b)

# #drop random quadratic
# model2 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus4_brainstem"))
# summary(model2)
# anova(model1,model2)
# #non-sig (marg), drop rand quadratic
# 
# model2b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus4_brainstem"))
# summary(model2b)
# anova(model1b,model2b)
# #non-sig (marg), drop rand quadratic
# 
# 
# #drop random slopes
# model3 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 | sub), data = subset(clpsc,roi=="zfclus4_brainstem"))
# summary(model3)
# anova(model2,model3)
# #non-sig, drop random slopes
# 
# model3b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 | sub), data = subset(clpsc,roi=="zfclus4_brainstem"))
# summary(model3b)
# anova(model2b,model3b)
# #non-sig, drop random slopes
# 
# 
# #remove fixed quadratic
# model4 <- lmer(psc ~ runc3*group*cond + (1 | sub), data = subset(clpsc,roi=="zfclus4_brainstem"))
# summary(model4)
# anova(model3,model4)
# #non-sig, drop fixed quadratic


model4b <- lmer(psc ~ runc3*group + (1 | sub), data = subset(clpsc,roi=="zfclus4_brainstem"))
summary(model4b)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: psc ~ runc3 * group + (1 | sub)
##    Data: subset(clpsc, roi == "zfclus4_brainstem")
## 
## REML criterion at convergence: -96.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9765 -0.5503 -0.0773  0.4900  3.9636 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  sub      (Intercept) 0.01396  0.1182  
##  Residual             0.02998  0.1732  
## Number of obs: 240, groups:  sub, 24
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)    0.156681   0.026583  22.000000   5.894 6.24e-06 ***
## runc3         -0.048754   0.007903 214.000000  -6.169 3.39e-09 ***
## group1         0.071485   0.026583  22.000000   2.689   0.0134 *  
## runc3:group1  -0.020254   0.007903 214.000000  -2.563   0.0111 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) runc3 group1
## runc3       0.000              
## group1      0.000  0.000       
## runc3:grop1 0.000  0.000 0.000
# anova(model3b,model4b)
# #non-sig, drop fixed quadratic

final_model=model4b

#plot by groups and day
g1 <- ggplot(data = subset(clpsc,roi=="zfclus4_brainstem"), aes(x = runc3, y = psc, shape = group, color=group))
g2 <- g1 + stat_summary(fun.data=mean_se, geom="pointrange", size=.5) + scale_colour_manual(values=c("#6666FF","#666666"))
g3 <- g2 + stat_summary(aes(y = fitted(final_model), linetype = group,color=group), fun.y = mean, geom="line")
g4 <- g3 + theme_bw(base_size=18) 
g5 <- g4 + scale_y_continuous(name = "Percent Signal Change") + scale_x_continuous(label=c("1","2","3","4","5"), name = "Session")
print(g5)

ggsave(filename="zfclus4_brainstem.png")
## Saving 7 x 5 in image

zfclus4_parahipp (LPHIPP)

# #full
# model1 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus4_parahipp"))
# summary(model1)

model1b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus4_parahipp"))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.511083 (tol =
## 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
summary(model1b)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: 
## psc ~ runc3 * group + I(runc3^2) * group + (1 + runc3 + I(runc3^2) |  
##     sub)
##    Data: subset(clpsc, roi == "zfclus4_parahipp")
## 
## REML criterion at convergence: -84.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.91306 -0.57881 -0.06313  0.60290  2.54873 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr       
##  sub      (Intercept) 0.0257529 0.16048             
##           runc3       0.0009823 0.03134   1.00      
##           I(runc3^2)  0.0018622 0.04315  -0.64 -0.64
##  Residual             0.0259395 0.16106             
## Number of obs: 240, groups:  sub, 24
## 
## Fixed effects:
##                    Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)        0.107183   0.036545 22.088000   2.933  0.00768 ** 
## runc3             -0.032220   0.009745 28.268000  -3.306  0.00258 ** 
## group1            -0.012812   0.036545 22.088000  -0.351  0.72921    
## I(runc3^2)        -0.006144   0.010779 22.108000  -0.570  0.57443    
## runc3:group1      -0.046757   0.009745 28.268000  -4.798 4.72e-05 ***
## group1:I(runc3^2)  0.023424   0.010779 22.108000   2.173  0.04077 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) runc3  group1 I(3^2) rnc3:1
## runc3        0.588                            
## group1       0.000  0.000                     
## I(runc3^2)  -0.668 -0.346  0.000              
## runc3:grop1  0.000  0.000  0.588  0.000       
## grp1:I(3^2)  0.000  0.000 -0.668  0.000 -0.346
## convergence code: 0
## Model failed to converge with max|grad| = 0.511083 (tol = 0.002, component 1)
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
# #drop random quadratic
# model2 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus4_parahipp"))
# summary(model2)

# anova(model1,model2)
# #sig, keep rand quadratic

model2b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus4_parahipp"))
summary(model2b)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: psc ~ runc3 * group + I(runc3^2) * group + (1 + runc3 | sub)
##    Data: subset(clpsc, roi == "zfclus4_parahipp")
## 
## REML criterion at convergence: -66.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.09581 -0.57346 -0.02772  0.59224  3.07387 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr
##  sub      (Intercept) 0.0148992 0.12206      
##           runc3       0.0006525 0.02554  1.00
##  Residual             0.0318926 0.17859      
## Number of obs: 240, groups:  sub, 24
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)         0.107183   0.030717  34.250000   3.489 0.001352 ** 
## runc3              -0.032220   0.009676  36.040000  -3.330 0.002014 ** 
## group1             -0.012812   0.030717  34.250000  -0.417 0.679197    
## I(runc3^2)         -0.006144   0.006889 212.000000  -0.892 0.373473    
## runc3:group1       -0.046757   0.009676  36.040000  -4.832  2.5e-05 ***
## group1:I(runc3^2)   0.023424   0.006889 212.000000   3.400 0.000805 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) runc3  group1 I(3^2) rnc3:1
## runc3        0.437                            
## group1       0.000  0.000                     
## I(runc3^2)  -0.449  0.000  0.000              
## runc3:grop1  0.000  0.000  0.437  0.000       
## grp1:I(3^2)  0.000  0.000 -0.449  0.000  0.000
anova(model1b,model2b)
## refitting model(s) with ML (instead of REML)
## Data: subset(clpsc, roi == "zfclus4_parahipp")
## Models:
## ..1: psc ~ runc3 * group + I(runc3^2) * group + (1 + runc3 | sub)
## object: psc ~ runc3 * group + I(runc3^2) * group + (1 + runc3 + I(runc3^2) | 
## object:     sub)
##        Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)    
## ..1    10 -88.603 -53.797 54.302  -108.60                             
## object 13 -99.779 -54.531 62.890  -125.78 17.176      3  0.0006502 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#sig, keep rand quadratic (but that model doesn't converge)

final_model=model2b

#plot by groups and day
g1 <- ggplot(data = subset(clpsc,roi=="zfclus4_parahipp"), aes(x = runc3, y = psc, shape = group, color=group))
g2 <- g1 + stat_summary(fun.data=mean_se, geom="pointrange", size=.5) + scale_colour_manual(values=c("#6666FF","#666666"))
g3 <- g2 + stat_summary(aes(y = fitted(final_model), linetype = group,color=group), fun.y = mean, geom="line")
g4 <- g3 + theme_bw(base_size=18) 
g5 <- g4 + scale_y_continuous(name = "Percent Signal Change") + scale_x_continuous(label=c("1","2","3","4","5"), name = "Session")
print(g5)

ggsave(filename="zfclus4_parahipp.png")
## Saving 7 x 5 in image

zfclus5 (LILOC)

# #full
# model1 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus5"))
# summary(model1)
# 
model1b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus5"))
summary(model1b)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: 
## psc ~ runc3 * group + I(runc3^2) * group + (1 + runc3 + I(runc3^2) |  
##     sub)
##    Data: subset(clpsc, roi == "zfclus5")
## 
## REML criterion at convergence: -129.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.55788 -0.50140 -0.02556  0.52108  3.03856 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr       
##  sub      (Intercept) 0.053545 0.23140             
##           runc3       0.000796 0.02821  -0.12      
##           I(runc3^2)  0.002243 0.04736  -0.48 -0.15
##  Residual             0.017461 0.13214             
## Number of obs: 240, groups:  sub, 24
## 
## Fixed effects:
##                    Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)        0.527696   0.049069 22.000000  10.754 3.16e-10 ***
## runc3             -0.050295   0.008339 22.000000  -6.031 4.53e-06 ***
## group1             0.047212   0.049069 22.000000   0.962 0.346424    
## I(runc3^2)         0.014331   0.010930 22.000000   1.311 0.203303    
## runc3:group1      -0.031926   0.008339 22.000000  -3.828 0.000916 ***
## group1:I(runc3^2)  0.005493   0.010930 22.000000   0.503 0.620229    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) runc3  group1 I(3^2) rnc3:1
## runc3       -0.083                            
## group1       0.000  0.000                     
## I(runc3^2)  -0.504 -0.093  0.000              
## runc3:grop1  0.000  0.000 -0.083  0.000       
## grp1:I(3^2)  0.000  0.000 -0.504  0.000 -0.093
# 
# 
# #drop random quadratic
# model2 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus5"))
# summary(model2)
# anova(model1,model2)
# #sig, keep rand quadratic

# model2b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus5"))
# summary(model2b)
# anova(model1b,model2b)
# #sig, keep rand quadratic

final_model=model1b

#plot by groups and day
g1 <- ggplot(data = subset(clpsc,roi=="zfclus5"), aes(x = runc3, y = psc, shape = group, color=group))
g2 <- g1 + stat_summary(fun.data=mean_se, geom="pointrange", size=.5) + scale_colour_manual(values=c("#6666FF","#666666"))
g3 <- g2 + stat_summary(aes(y = fitted(final_model), linetype = group,color=group), fun.y = mean, geom="line")
g4 <- g3 + theme_bw(base_size=18) 
g5 <- g4 + scale_y_continuous(name = "Percent Signal Change") + scale_x_continuous(label=c("1","2","3","4","5"), name = "Session")
print(g5)

ggsave(filename="zfclus5.png")
## Saving 7 x 5 in image

zfclus6 (RILOC)

# #full
# model1 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus6"))
# summary(model1)

model1b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus6"))
summary(model1b)
## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: 
## psc ~ runc3 * group + I(runc3^2) * group + (1 + runc3 + I(runc3^2) |  
##     sub)
##    Data: subset(clpsc, roi == "zfclus6")
## 
## REML criterion at convergence: -82
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3977 -0.5061 -0.0496  0.5509  2.6569 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr       
##  sub      (Intercept) 0.033945 0.18424             
##           runc3       0.001953 0.04419   0.24      
##           I(runc3^2)  0.001007 0.03173  -0.34 -0.54
##  Residual             0.023633 0.15373             
## Number of obs: 240, groups:  sub, 24
## 
## Fixed effects:
##                    Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)        0.537071   0.040663 22.000000  13.208 6.17e-12 ***
## runc3             -0.066311   0.011428 22.000000  -5.803 7.74e-06 ***
## group1             0.050298   0.040663 22.000000   1.237   0.2292    
## I(runc3^2)         0.013941   0.008782 22.000000   1.587   0.1267    
## runc3:group1      -0.031577   0.011428 22.000000  -2.763   0.0113 *  
## group1:I(runc3^2)  0.005463   0.008782 22.000000   0.622   0.5403    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) runc3  group1 I(3^2) rnc3:1
## runc3        0.172                            
## group1       0.000  0.000                     
## I(runc3^2)  -0.430 -0.312  0.000              
## runc3:grop1  0.000  0.000  0.172  0.000       
## grp1:I(3^2)  0.000  0.000 -0.430  0.000 -0.312
# #drop random quadratic
# model2 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus6"))
# summary(model2)
# anova(model1,model2)
# #sig, keep rand quadratic
# 
# model2b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus6"))
# summary(model2b)
# anova(model1b,model2b)
# #sig, keep rand quadratic


final_model=model1b

#plot by groups and day
g1 <- ggplot(data = subset(clpsc,roi=="zfclus6"), aes(x = runc3, y = psc, shape = group, color=group))
g2 <- g1 + stat_summary(fun.data=mean_se, geom="pointrange", size=.5) + scale_colour_manual(values=c("#6666FF","#666666"))
g3 <- g2 + stat_summary(aes(y = fitted(final_model), linetype = group,color=group), fun.y = mean, geom="line")
g4 <- g3 + theme_bw(base_size=18) 
g5 <- g4 + scale_y_continuous(name = "Percent Signal Change") + scale_x_continuous(label=c("1","2","3","4","5"), name = "Session")
print(g5)

ggsave(filename="zfclus6.png")
## Saving 7 x 5 in image

plot all together

clpsc$roi<-factor(clpsc$roi,levels=c("zfclus1_lifg","zfclus1_lins","zfclus2","zfclus3","zfclus4_brainstem","zfclus4_parahipp","zfclus5","zfclus6"),labels=c("LIFG", "LINS","LSLOC","RSLOC","BSTEM","LPHIPP","LILOC","RILOC"))
clpsc$roi<-factor(clpsc$roi,levels=rev(levels(clpsc$roi)))

#facet grid with freq only
g1 <- ggplot(data = subset(clpsc,roi != "NA" & cond=="freq"), aes(x = as.numeric(run), y = psc, colour=group))
g2 <- g1 + stat_summary(aes(y = psc,linetype=group,color=group), size=.5, fun.y = mean, geom="line") + stat_summary(fun.data=mean_se, geom="pointrange", size=.5) + scale_colour_manual(values=c("#6666FF","#666666")) + scale_colour_manual(values=c("#6666FF","#666666")) + facet_grid(~roi)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
g3 <- g2 + theme_bw(base_size=12)  
g4 <- g3 + scale_x_continuous(name = "Session") +scale_y_continuous(name = "Percent Signal Change") 
print(g4)

ggsave(filename="facet_grid_pscs_ROIs.png",width=7,units=c("in"),dpi=300)
## Saving 7 x 5 in image
#facet grid with freq & infreq 
g1 <- ggplot(data = subset(clpsc,roi != "NA"), aes(x = as.numeric(run), y = psc, colour=group))
g2 <- g1 + stat_summary(aes(y = psc,linetype=group,color=group), size=.5, fun.y = mean, geom="line") + stat_summary(fun.data=mean_se, geom="pointrange", size=.5) + scale_colour_manual(values=c("#6666FF","#666666")) + scale_colour_manual(values=c("#6666FF","#666666")) + facet_grid(~roi*cond)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
g3 <- g2 + theme_bw(base_size=12)  
g4 <- g3 + scale_x_continuous(name = "Session") +scale_y_continuous(name = "Percent Signal Change") 
print(g4)

ggsave(filename="facet_grid_pscs_withINF_ROIs.png",width=10,units=c("in"),dpi=300)
## Saving 10 x 5 in image